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1. Motivation and Objective 

Recent development using compact difference schemes to solve the Navier-Stokes 
equations show spectral-like accuracy 1 - 2 . In this paper, we report further study 
of the numerical characteristics of various combinations of the Runge-Kutta (RK) 
methods and compact difference schemes to calculate the unsteady Euler equations. 
Conventionally, the accuracy of finite difference schemes is assessed based on the 
evaluations of dissipative error. The objectives are reducing the numerical damp- 
ing and, at the same time, preserving numerical stability. While this approach has 
tremendous success solving steady flows, numerical characteristics of unsteady cal- 
culations remain largely unclear. For unsteady flows, in addition to the dissipative 
errors, phase velocity and harmonic content of the numerical results are of concern. 
As a result of the discretization procedure, the simulated unsteady flow motions 
actually propagate in a dispersive numerical medium. Consequently, the dispersion 
characteristics of the numerical schemes which relate the phase velocity and wave 
number may greatly impact the numerical accuracy. The objective of the present 
paper is to assess the numerical accuracy of the simulated results. To this end, 
the Fourier analysis is performed to provide the dispersive correlations of various 
numerical schemes. 

First, a detailed investigation of the existing RK methods is carried out. A 
generalized form of an N-step RK method is derived. With this generalized form, 
the criteria are derived for the three and four-step RK methods to be third and 
fourth-order time accurate for the non-linear equations, e.g., flow equations. These 
criteria are then applied to commonly used RK methods such as Jameson s 3- 
step and 4-step 3 ’ 4 schemes and Wray’s algorithm 5 to identify the accuracy of the 
methods. For the spatial discretization, compact difference schemes are presented. 
The schemes are formulated in the operator-type 6 to render themselves suitable for 
the Fourier analyses. The results of the analyses provide CFL limits, the numerical 
dispersion relations, and the artificial damping required for stable and time-accurate 
solutions. 

Finally, the performance of the numerical methods is demonstrated by numeri- 
cal examples. The first case is a quasi-one-dimensional calculation of the acoustic 
admittance in a converging nozzle. The CFD results are compared with Tsien’s 
analytical solution 7 ; the harmonic content of this flow field is limited to one fre- 
quency mode. All numerical schemes of concern provide accurate solutions. The 
second case is a one-dimensional simulation of a shocked sound wave. The harmonic 




127 

PA'iE .55LA5M 


« a at 





S.-T. Yu 


content is complex and distinct differences between various schemes are observed. 
The results are also compared with the analytical solution provided by Morse and 
Ingard 8 . In the one-dimensional cases, details of the numerical methods in setting 
up the initial conditions and the perturbation on the computational boundary are 
described. 

The third case is a two-dimensional simulation of a Lamb vortex 9 in an uniform 
flow. This calculation provides a realistic assessment of various finite difference 
schemes in terms of the conservation of the vortex strength and the harmonic content 
after travelling a substantial distance. The numerical implementation of Giles’ non- 
reflective equations 17 coupled with the characteristic equations as the boundary 
condition is discussed in detail. Finally, the single vortex calculation is extended 
to simulate vortex pairing 10 . For the distance between two vortices less than a 
threshold value, numerical results show crisp resolution of the vortex merging. 


2. Work Accomplished 
2.1 Numerical Method 

The Euler equations in Cartesian coordinates can be cast into a vector form: 


5Q f 9Ei 

dt dx{ 

1=1 


= 0 , 


(i) 


where Q is the unknown vector and E* is the inviscid flux in the X{ direction. 
The Runge-Kutta algorithm is applied as the temporal discretization and the sec- 
ond, fourth, and sixth-order compact difference schemes are applied to the spatial 
discretization. 


2.2.1 The Runge-Kutta Method 

The use of the Runge-Kutta methods for flow equations stems from the applica- 
tion of the methods to solve ordinary differential equations (ODEs). An ODE has 
one independent variable and its solution is obtained by integrating the equation 
from its initial condition. When one applies the Runge-Kutta method to the flow 
equations, time is treated as the independent variable as is in an ODE, and the 
convective terms are taken as the inhomogeneous part of the equations, such as 


w m R(Q) • <2) 

Notice that the boldface symbols which represent vectors have been temporarily 
dropped for typographic convenience. In addition, all results in the following dis- 
cussion are valid for both scalar and vector equations. 

The Runge-Kutta methods have algorithms of the form 


Q n+ 1 = Q n + At R(Q n , At), 


( 3 ) 
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where the increment function R(Q n ,At) is a suitable chosen approximation to the 
inhomogeneous part of the equation, that is, R(Q)- In general, t e c cu a ion+0 
the increment function R is subdivided into N steps on the interval t <t _t 
And the final increment function R is a weighted average of the inhomogeneous 
terms evaluated at the different steps on the interval t n <t <t n , that is 

Q l = Q n + At{a n R n ), 

Q 2 = Q n + At(a 2 \R n + £*22 R 1 ), 

Q z = Q n + At(a 3 iR n + £* 32 #* + £*33^ 2 )» (4) 


Q n+1 = Q n + At(otNiR n + OtmR 1 + • • ■ + OtffNR N 1 ), 

where the superscript n, 1, 2, .... and n+1 denote the time steps on the time interval 
t n < t\ < t 2 < ... < ttf < < n+1 , and ctij is the weighting factor for the step i and 

term j. There are i weighting coefficients to be determined and an infinite 

number of coefficient sets can be chosen. However, certain criteria must be met for 
the algorithm to retain high-order accuracy. 

In what follows, the criteria of the coefficient set of a 3-step Runge-Kutta method 
to be third-order accurate is given. To proceed, we follow the conventional approach 
and expand all inhomogeneous terms R * in Eqn. (4) to a Taylor’s series about R 
and drop all terms in which the exponent of At is greater than 3. The result 
is compared with its analytical counterpart by equating terms m like powers of 
At. The result is tabulated in Table 1. For the convenience of the discussion, the 
following simplification of symbols is activated: R denotes R(Q n ), Q denotes Q n , 
R' denotes {dR/dQ) n , and R" denotes ( d 2 R/dQ 2 ) n . In addition to the equality of 
the coefficients of all the powers of At, we also want the equality of the coefficients 


tug v Vill V1VU VW V* X * 

of the functions of R, R', and R" . As a result, we find the 
for the 3-step RK methods to be third-order accurate as, 

criteria of the coefficients 

£*31 + £*32 + £*33 = 

(5a) 

£*11 £*32 + £*33 (£*21 + £*22) = 2’ 

(56) 

, \2 1 
C*j j £*32 + (<*21 + <*22) <*33 — j? 

(5c) 

1 

£*11<*22£*33 — g- 

( 5<f) 


Equations (5a) and (5b) are the criteria of first and second order accuracy, re- 
spectively. The remaining equations are of third order term. Since four equations 
contain six unknowns, the system is underdetermined, and two of the coefficients 
may be chosen arbitrarily. The obvious choice is to let the two coefficients be null 
to reduce intermediate storage and numerical operations. According to Eqn. (5d), 
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none of ecu where i = 1, 2, 3, could be zero, and one can set the two of the three re- 
maining coefficients to be zero. Therefore, at least one of the intermediate steps has 
two non-zero coefficients. Consequently, one needs to store two steps of intermediate 
solutions for the 3-step, third-order RK methods. 

A 3-step Runge-Kutta method proposed by Jameson et al. 3 to solve flow equations 
is 


Q 1 = Q n + A tR n , 

= + + ( 6 ) 
<?" +1 = Q” + y(R" + R 3 ). 

It can be shown that the weighting coefficients of the 3-step method satisfy only 
Eqns. (5a) and (5b). And the method is second-order accurate in time. 

Wray 5 proposed another 3-step method, 

«■=(?’> + 

Q 3 = Q l + At(~R" - ^R 1 ), ( 7 ) 

<? n+1 = e ! + A((jfi"- ^R 3 ). 

This formula may be manipulated to fit the generalized form as proposed in Eqn. 
(4), and we obtain, 

<?'=< 3" + AI(lfl"), 

<? ! = <?" + A((jfl" + ^fl 1 ). (8) 

4 4 

Wray’s coefficients match all the equations in Eqn. (5) and therefore the scheme is 
third-order accurate. In this formula, <*32 is set to zero and two sets of solutions 
are needed in the second and the third steps. The calculation can be carried out by 
either the vectorized algorithm proposed by Wray, or straightforward calculation 
according to Eqns. (7) and (8). 

A similar procedure can be applied to the 4-step Runge-Kutta methods, and the 
criteria for the scheme to be fourth-order accurate are: 


Q41 + &42 + 0:43 + <*44 = 1, (9a) 

&lia42 + (0:21 + <*22)0:43 + (<*31 + <*32 + <*33)044 = 2’ (9&) 
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A _ 1 

< 21 x 0*42 + (<*21 + <* 22) 2 <*43 + (<*31 + <*32 + <*33 ) <*44 — g i 

(9c) 

<*11 <*22 <*43 + [<*11 <*32 + (<*21 + <* 22 )<* 33 ] <*44 = gf 

( 9 d ) 

\3 1 

<* 11«42 + (<*21 + <* 22 ) 3 <*43 + (<*31 + <*32 + <*33 ) <*44 - 

(9e) 

-<* n <* 22<*43 + (<*21 +<* 22 )<* 11 <* 22<*43 
+ ^ [(<*21 + <*22 ) 2 <*33 + <*11 <* 32 ] <*44 


+ (<*31 +<*32 + <* 33 ) [<*11 <*32 + (<*21 +<* 22 )<* 33 ] <*44 = gi 

(9/) 

1 

<*11 <*22 <*33 <*44 — 24 ' 

(99) 


Equations (9a) and (9b) are for first and second-order accuracy, respectively. Equa- 
tions (9c) and (9d) are for third-order accuracy. The remaining equations are for 
the fourth-order terms. Here, seven equations contain ten unknowns, and three of 
the coefficients may be chosen arbitrarily. 

A 4-step RK method attributed to Kutta 11 for solving ODEs was adopted by 
Jameson et al 3 to solve the flow equations. The algorithm can be expressed as, 

q 1 = cr + y r "’ 

<? 2 =Q" + Y*‘' (io) 

Q 3 = Q n + AtR 2 , 

gn+l + 2 R l + 2 R 2 + R*)- 

The coefficients satisfy Eqn. (9) and the algorithm is fourth-order accurate. How- 
ever, this method requires all four intermediate solutions in the final step. As a 
result, the use of this scheme for large-scale calculations is undesirable. 

Later on, Jameson and Baker 4 proposed another 4-step algorithm, 

Q l = Q n + y R n , 

C 2 = <3" + T«'- (ii) 

«’ = <3" + Y« 2 . 

Q n +* - Q n + AtR 3 . 

This scheme was proposed to calculate the steady state solutions and the transient 
solutions were not of concern. For that purpose, the algorithm is convenient to 
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program and no intermediate solution needs to be stored. For the present investi- 
gation, however, the weighting coefficients of this method satisfy only part of Eqn. 
(9) and the algorithm is second-order accurate. Nevertheless, this formulation is 
favorable compared to a 2-step, second-order RK method because part of the third 
and fourth-order terms are satisfied, namely, Eqns. (9d) and (9g). Consequently, a 
larger marching step, i.e., a larger CFL number, could be used. 

2.2.2 Compact Difference Schemes 

The remaining task of discretizing the flow equations is the spatial differencing 
of the inviscid fluxes. An effective manner for generating a high-order, central 
difference scheme is the compact difference method. The scheme is obtained by 
using only three and five points to achieve fourth-order and sixth-order accuracy 
in space, respectively. The gain in the accuracy is not based on the involvement of 
more points as in the conventional approach, but on implicitly solving the derivatives 
simultaneously at all locations. According to Hermite’s generalization of a Taylor’s 
series, 12 one can get 


Ui_ x + 4 u'i + = — (u i+1 - Mi- 1) + 0(Ax 4 ), (12) 

Mi_i + + u' i+1 = (m*+ 2 + 28u i+ i - 28«i-i - Mi- 2 ) + 0(Ax 6 ), (13) 

where the superscript ' represents the spatial derivatives. Equation (12) is the 
fourth-order method and Eqn. (13) is the sixth-order one. When the fourth- 
order method is used in the interior nodes, a third-order biased implicit scheme 13 
is adopted for grid nodes at the computational boundary, such as 

2u 1 + 4*4 = — ^-(-5i/i -1- 4 U 2 + u 3 ) + 0(Ax 3 ), 

1 ( 14 ) 
^M moz -(- 4u maz _j — (5u mai 4u maz _i — u mai _2) + 0(Ax 3 ). 

When the sixth-order method is used in the interior nodes, the fourth-order scheme, 
Eqn. (12), is used at locations one grid node away from the boundary and the 
third-order biased scheme is used at the boundary. The application of the implicit 
compact difference schemes with the appropriate boundary conditions involves the 
inversion of a scalar tridiagonal matrix. The inversion of the matrix incurs little 
penalty in terms of CPU time. 

2.3 Fourier Analysis 

By definition, any function, u(x,t), which is continuous, periodic, and square 
summable can be expressed in a Fourier series expansion at a constant time, 


OC 

u(x,t)= ^ u(k,t)e t2nkx / L (15) 

k= — 00 
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where L is the period of the function u(x,f), k is the wave number, and i is -/-l- 
The Fourier coefficient is defined as, 

u(k,t) = 7 f ' u(x,t)e~ i2 * kx/L dx. (16) 

L J-L/2 

In the Fourier analysis of a finite difference scheme, functions are defined at discrete 
points. The discrete Fourier series and its coefficients are defined analogous to their 
continuous counterparts, 


K - 1 

u] = u n (k)e“* kj/K , j = o, ±1, ±2, • ■ • , ± 00 , 

fc=o (17) 

u"(h) = ± £ u?e- i2 *W K , k = 0, 1, • • • , K - 1. 
j= i 

Here, the length of the computational domain L is decomposed into K grid nodes 
(£, = K Ax). The superscript n denotes the time step and the subscript j is the 
spatial location index. Similar to the continuous function, the algebraic system in 
terms of the function exp(i2nkj/K) is periodic over the computational domain L 
(or K ) and orthonormal, such as 

gilnkj/K _ gi2*k(j+K)/K 

1 K ~ l (18) 

K r-i 

where fikik 2 I s the Kronecker delta. Therefore, the establishment of the discrete 
Fourier series and coefficients is self-sufficient, and is not an approximation of its 
continuous counterpart. 

As shown in Eqns. (17) and (18), the harmonic content of the discretized equation 
is limited to the number of grid nodes used in the computational domain. A discrete 
solution u] at a location (j) and time (n) is a Unear combination of K wave modes. 
The Fourier analysis is performed by substituting each wave mode of the discrete 
Fourier expansion, Eqn. (17), into the discretized flow equations to calculate the 
amplification factor, g(k), which is defined as 


9( k ) = 


u n+1 (fc) 

u n (k) 


(19) 


The procedure is repeated for all wave modes (fc = 0,1, ••• ,K - I) and the full 
spectrum of the amplification factor is obtained. In this process, we map the func- 
tion u in terms of spatial variable x on the interval [-L/ 2 ,L/ 2 ] to the wave number 
space on [— 7 r, 7 r] assuming that the analysis is local for an infinite and periodic do- 
main. Therefore, the solution of the amplification factor on the interval [-7T, 0] is 
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the complex conjugate of that on [0,7r], For this reason, the results of our Fourier 
analyses are presented on the interval [0,7r]. 

In the present investigation, one-dimensional equations are considered for the 
Fourier analysis. In addition, we can perform a similarity transformation to trans- 
form the one-dimensional Euler equations to their characteristic form, i.e., three 
decoupled scalar equations. Consequently, a scalar, advective equation on a peri- 
odic domain is adopted as the model equation in our analysis, 


du du 
~dt + *di 


= 0 , 


( 20 ) 


where the phase velocity (A) is equivalent to the eigenvalues of the Euler equations, 
namely u - c, u + c, and u where u is velocity and c is the speed of sound. The 
phase speed (A) is treated as a parameter in the Fourier analysis to avoid the 
Fourier convolution, therefore the analysis is linear. For the unsteady calculations, 
the requirement of the time resolution of the flow field restricts the time marching 
step. In other words, the variations of flow properties between time steps are small. 
Thus, linear analysis is a viable tool. 


In what follows, the procedure to obtain the amplification factor of the model 
equation discretized by Runge-Kutta methods and compact differences is illustrated. 
First, the generalized forms of the amplification factor for the 3 and 4-step Runge- 
Kutta methods are derived. These representations of the amplification factors are 
independent of the spatial discretization schemes. From the equations of Wray’s 
3-step scheme, Eqn.(7), we have 


9 2 — 1 + + Y2^g l , (21) 

9 = 1 + \ Z + \ Z 9 2 ' 


where 


g 1 Q n = Q 1 , 

g 2 Q n = Q 2 , ( 22 ) 

gQ n = Q n+1 . 

The variable Z represents the spatial discretization applied to the convective term ( 
-Xdu/dx ). By substituting the amplification factors of the intermediate steps, g 1 
and g 2 , into the last step of Wray’s algorithm, we obtain the amplification factor, 
g , for the 3-step scheme, 


S = l + Z + iz 2 + iz 3 . 


(23) 
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It is interesting to note that Eqn. (23) can be directly derived from the Taylor s 
series expansion by adopting the invariance property of the time and spatial deriva- 
tives of the model equation, i.e., Z = —A du/dx = dufdt. This is valid because the 
coefficients satisfy Eqn. (5), which is deduced from the Taylor’s series expansion up 
to the third-order term. On the other hand, the amplification factor of the 3-step 
scheme proposed by Jameson et al., Eqn. (6), can be derived as 

g = l + Z+lz 2 + ]z\ (24) 

* 2 4 

It is obvious that the scheme is not third-order accurate. 

A similar analysis can be applied to the 4-step methods. Identical forms of the 
amplification factors of both 4-step methods of concern (Eqns. (10) and (11)) are 
obtained, such as 

g=l + Z + \z'‘ + \z 3 + ^Z > . (25) 

Unlike the case of the 3-step schemes, the effect of the order of accuracy of these 
two 4-step schemes does not appear in the expression of the amplification factor. 
This is because the amplification factor is derived based on the linearized equation. 
Only by using Eqns. (5) and (9), which take into account nonlinear terms, can one 
justify the order of the accuracy of the RK schemes. 

The remaining task is to derive the explicit form of the spatial discretization 
operation, Z, of compact difference schemes. The fourth-order compact difference 
method, Eqn. (12), can be cast into the operator-type by defining 

6 2 Ui = Ui+\ — 2ui + Uf-i, (26) 

where Ui could be any flow property of interest at grid point i. As a result, the 
fourth-order method can be rewritten as, 

)+0(Al4) - (27) 

This equation allows us to express ( dujdx)i in an explicit form, 

(£),4 + £H"^) +0(Ai4) - (28) 

To proceed, we substitute this explicit, discretized form, Eqn. (32), into the model 
equation, and we obtain 


Z (4) = 


6Fsin(k)i 
4 + 2cos(k)' 


(29) 


135 


S.-T. Yu 


where F is the CFL number which is defined as F = XAt/Ax , and k is the normal- 
ized wave number (k — 2nk/K) . 

It is interesting to note that if the solution reaches a steady state solution, i.e., the 
time derivative term is zero, the operator (1 + 6 2 /6)~ l in the discretized equation 
becomes futile and the spatial discretization is represented by («i+i — ti,_i)/2Ax, 
which is only second-order accurate. However, the steady state solution of the one- 
dimensional wave equation is a constant and the spatial accuracy is meaningless. 
On the other hand, the accuracy of multi-dimensional calculations is more complex. 
For example, consider a two-dimensional version of the model equation discretized 
by the fourth-order compact difference method, and we have 


du x ( 

w + M 


, 

1+ i 


-1 


\ 2Ax 




(30) 

Again, the operators (1 + 5*/6) and (1 + 6 2 /6) are scalar tridiagonal matrices with 
the dimensions IL x IL and JL x JL, respectively. IL and JL are the numbers 
of the grid nodes in the x and y directions of the computational domain. When 
IL = JL, two operators are identical. We then multiply Eqn. (30) by the operator 
(1 + 6 2 f 6) and obtain the steady state equation as 


\ u i+l,j , x u i,j+l u i,j - 1 __ n 

Al 2Ai +A » 2A ~y 

Therefore, the steady state solution is only second-order accurate. 


( 31 ) 


Similarly, the spatial discretization of the sixth-order compact difference scheme 
can be represented in the operator-type such as, 


du 

dx 



1 

60Ax 


(it,: + 2 + 28tii + i - 28ui-i - m- 2 ) + 0(Ax 6 ). 


(32) 


And we obtain, 


= 


F[4sin(fc)cos(fc) + 56sin(fc)]i 
12[2cos(fc) + 3] 


(33) 


Again, when solving the steady state solution with same number of grid nodes in 
each direction of the computational domain, the spatial discretization is represented 
by (ui+2+28ui + i~28ui-i -«j_2)/60Ax. Unlike the case of the fourth-order scheme, 
this representation does not match any conventional central difference scheme and 
it is at most second-order accurate. 


According to the above discussion, we can obtain the amplification factor g(k) for 
various combinations of the Runge-Kutta methods and compact difference schemes. 
The amplification factor g(k) is a complex number and can be expressed as g(k) = 
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exp[iui(k)], where u) is the normalized frequency and is defined as u = 27ru;A tjr, 
u is the frequency, r is the the time period of function u and the phase speed 
(A) is equal to L/r. Here, the dissipative and dispersive artifacts of the numerical 

schemes can be assessed: 

1. Dissipation. The normalized frequency u is a complex number (d; = a + ip) and 
its imaginary part represents the magnitude of the amplification factor, i.e., 

9(k) = = e-'V” (34) 

|s(£)| = e~ 0 

The magnitude of the amplification factor is the artificial dissipation. When 
I0I > the scheme is unstable. For the calculations of unsteady flows, we want 
\ g \ to be less than unity but very close to it to ensure numerical stability with 
minimum artificial dissipation. In the following section, we plot \g\ against k to 
illustrate the artificial dissipation. 

2. Dispersion. According to Eqn. (34), the relation of Q(fc) and k represents the 
artificial dispersion. We plot & = a/F against k to show phase velocities. Notice 
that the model equation is dispersionless and the phase velocity is a constant, i.e., 
A. After being normalized by the CFL number, the exact solution is a straight 
line with 45° angle on the plot of & against k. 

Figure 1 shows the results of the Fourier analysis of the third-order Runge-Kutta 
(RK3) method combined with fourth (CD4), and sixth (CD6) order compact differ- 
ence schemes and the conventional second-order central difference scheme (CD2). 
The figures show the dissipative as well as dispersive effects at CFL numbers from 
0.4 to 1.4 with an increment of 0.2 between neighboring curves. Figure la shows 
the dissipation of the RK3-CD6 scheme. The method is unstable for CFL numbers 
greater than 0.8. As the order of spatial differencing decreases (compare Figs, la, 
lc, and le), the limit of the CFL number increases for stable calculation. 

Figure lb shows the dispersive effect of the RK3-CD6 scheme. For a CFL number 
of 0.4, the phase velocity is correct for wave numbers up to 1.8. The phase velocities 
are slower than they should be at large wave numbers. Increasing the CFL number 
makes phase velocities deviate from the 45° straight line at smaller wave numbers. 
Decreasing the CFL number merges the curves together to reach an asymptotic 
curve. However, the dispersive effect at high wave numbers does not improve. 
Comparing Figs, lb, Id, and If shows that the increase of the order of the spatial 
differencing reduces the numerical dispersion at high wave numbers. Specifically, a 
significant improvement is achieved by changing the spatial differencing from CD2 
to CD4, whereas only a limited gain is obtained by switching from CD4 to CD6. 

Figure 2 shows the results of the Fourier analyses of the fourth-order Runge-Kutta 
(RK4) method combined with various central difference schemes. The CFL numbers 
are the same as that in Fig. 1. Similar to the case of RK3, for the same CFL number 
and wave number, e.g., CFL = 1.4, u> = 1.5, higher-order spatial discretization 
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introduces more artificial damping (see Figs. 2a, 2c, and 2e) and therefore reduces 
the CFL number limit for stable calculation. Again, the dispersive error at high 
wave numbers decreases as the order of the spatial differencing increases (see Figs. 
2b, 2d, and 2f). 

Figures, lc and 2c can be compared to show the difference of the dissipation 
effects between the RK3 and RK4 methods. For the same CFL and wave numbers, 
the RK4 method introduces more artificial damping, and a larger CFL number 
could be used. On the other hand, Figs. Id and 2d show that an increase of the 
order of the time marching scheme does not improve the dispersive effect at high 
wave numbers. 

From the above discussion, it is clear that reliable solutions of the finite difference 
schemes are at low wave numbers. For example, for the RK4-CD6 method at CFL 
= 0.8 (see Figs. 2a and 2b) the solution with wave numbers less than 1 /3tt (6 grid 
nodes per wave) is fairly accurate. Numerical solutions with higher wave numbers 
(wave length less than 6 grid nodes) suffer significant dispersive and dissipative 
errors. On the other hand, for the conventional RK4-CD2 method at the same 
CFL, 12 to 16 grid nodes per wave are needed for an accurate solution. 

It is interesting to note that compact difference schemes have no dissipative effect 
at the highest wave number resolved by a given numerical grid, i.e., two grid nodes 
per wave (see Figs, la, lc, le, 2a, 2c, and 2e). Nevertheless, a significant dispersive 
error is introduced to these highest- wave-number waves and cause the even-odd de- 
coupling of the numerical solutions. Furthermore, applying the compact difference 
scheme twice to calculate the viscous terms of the Navier Stokes equations does not 
eliminate the erroneous oscillation, owing to the linearity of the operation. These 
high-wave-number waves continue oscillating with erroneous phase speeds through- 
out the course of computation and eventually destroy the solution. It is therefore 
appropriate to impose a small amount of high-order artificial damping to filter out 
these waves while at the same time keeping the resolution at low wave modes intact. 
Figure 3 shows the dissipation and dispersion effects of the RK4-CD6 method at 
CFL=0.8 with various amounts of sixth-order artificial damping, defined as 

A.D. = J [u i+3 + Ui- 3 - 6{u i+2 + Ui-i) + 15(ll i+ 3 + Ui - 3 ) - 20u<]. (35) 

O 

The range of T) is from 0.01 to 0.05 with an increment of 0.01 between the neigh- 
boring curves. Comparison between Fig. 3 and Figs. 2a and 2b shows that no 
additional damping at low wave numbers is introduced into the system by the im- 
posed artificial dissipation for tj < 0.03, whereas the undesirable high-wave-number 
waves are dissipated. 

2.4 Numerical Examples 

2.4.1 Acoustic Admittance of A Nozzle Flow 

The first case is a forced oscillatory quasi-one-dimensional flow in a converging 
nozzle. The governing equations are 
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dQ dE 
dt dx 


= H 


(36) 


where 

E= ^ pj + p j a, H = ^J> 3 § j , ( 37 ) 

p is density, p is pressure, and e is the total energy defined as e = p(C y T + \u 2 ). C v 
is the constant volume specific heat. The variable a is the cross sectional area and 
is prescribed as a function of x. The theoretical solution of the acoustic admittance 
of a choked nozzle was provided by Tsien 7 under the assumption that the velocity 
of the base flow is a linear function of axial location as shown in Fig. 4. The nozzle 
shape can be inversely derived according to Tsien’s assumption, and we have 


1 




+ 


7 ~ 1 
7 + 1 



2 (^- 1 ) 


3 


( 38 ) 


where 7 is the specific heat ratio and M is the Mach number which can be expressed 
as 



7-1 ( M 2 

2 \X*) 


( 39 ) 


The superscript * denotes the property at the nozzle throat. According to Tsien s 
derivation, the linearized quasi-one-dimensional equations can be manipulated to 
the following form under the isentropic condition, 


,d?P J, , iP \ dP (2 + 0) 

!(1 - z) l? ( + 1T7 i _ /J 207+Tj P _ 0 


dP 


(7 + !)(! — — (7 — 1 *+* ifi)P + (2 + 0W ~ 0 


where 


( 40 ) 

( 41 ) 


£- = P(z)e i0T , 

IV (42) 

- = U(z)e i0T . 
u 

u and p are the velocity and pressure of the base flow, 0 is the normalized frequency 
which is defined as 0 = w(l - z)/(c* - u ), and r is the non-dimensionalized time 
which is defined as r = c*t/x*. The independent variable z can be expressed 
in different forms due to the linearity between the base flow velocity u and axial 
location x, and we have 
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x 



u 


= ^ (43) 

( 7 + l)M 2 

2 + (7 - 1 )M 2 

With Eqn. (43), it is clear that P and U are functions of the Mach number (M) 
with prescribed frequency (/?). Equation (40) is a hypergeometric equation 14 that 
can be solved by a power series expansion. The converged solution does not exist 
in the supersonic region because the Mach number is greater than unity. U(z) can 
be easily solved with P(z) known as shown in Eqn. (41). Finally, the acoustic 
admittance function defined as A(z ) = U(z)/P(z ) can be obtained as a function of 
the Mach number. 

In what follows, the procedure of the CFD calculation to compare with Tsien’s 
solution is illustrated. First, the base flow field is obtained by solving the quasi- 
one-dimensional equations, Eqns. (36) and (37), using the RK4-CD2 method with 
the nozzle area ratio prescribed by Eqns. (38) and (39). The results are checked 
by the area Mach number relation 15 and the solution is accurate up to five decimal 
digits. The perturbation at the inlet is obtained by specifying sinusoidal pressure 
fluctuations in terms of magnitude and frequency. With the prescribed pressure 
and isentropic correlation, the temperature fluctuation is also determined. Numer- 
ically, these boundary conditions are enforced by defining a vector k = k(Q) at the 
upstream boundary, such as 


k = 



(44) 


where £1 and {2 are the specified values of p and T. To proceed, Equation (44) is 
linearized to become a function of AQ, such as 


k n+1 = k" + ^ AQ, (45) 

where k n+1 is equal to the specified pressure and temperature at the time step 
n + 1 and dk/dQ is a 3 x 3 matrix. In order to close the system, the null entry 
in the vector k may be filled by the out-running characteristic relation deduced 
from the flow equations. Numerically, the similarity transformation is applied to 
the discretized flow equations (see Eqn. (4)), and we get 


LM _1 (Q’ - Q n ) =LM- 1 At^a^.R fc - 1 , 

fc=i 


i = !,••' t N 


(46) 
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where i = 1, • • ■ , JV represents the N-step RK method. Here,- M _1 is the eigenvector 
matrix of Jacobian matrix A = dE/dQ, and L is a selection matrix with zeros and 
ones on the diagonal in such a fashion that the proper outrunning characteristics 
are selected. By combining the imposed conditions, Eqn. (45), with the outrunning 
characteristic relations, Eqn. (46), we form the complete equation at the boundary 

point as, 

(lm- 1 + H) (Q* - Q”) = LM-'AlX^R*- 1 +k" +I -k", (4?) 

i = 1 ,- ,N 

For the supersonic out-flow condition, Eqn. (46) is used with the selection matrix L 
equal to an identity matrix. In both cases, the out-running characteristic equations 
are solved with one-sided difference as shown in Eqn. (14). In other words, the 
characteristic boundary conditions are always discretized by an upwinding scheme 
which is physically sound and the numerical stability is enhanced. These boundary 
conditions are applied at each of the Runge-Kutta stages. 

The acoustic admittance is a complex number and can be written as A = \A\e l9 . 
In the present paper, a small pressure perturbation of 1.1% ( p ' = 0.01 lp) is imposed 
at the nozzle inlet. The length of the converging part of the nozzle is 0.9 L* (see 
Fig. 4) and the inlet Mach number is about 0.09. The frequency of the perturbation 
is set at ft = 6, which corresponds to about 2000 Hz. 

Figure 5 shows the comparisons between the CFD results of the RK4-CD6 method 
and the theoretical solution of the acoustic admittance in terms of the magnitude 
\A\ and the phase angle 0 in the subsonic region of the nozzle. Both the magnitude 
and the phase angle of the acoustic admittance decrease as the flow speeds up. As 
shown in the figure, perfect agreement is obtained for the comparison of |A|, while 
the predicted phase angle is slightly off due to the resolution of the numerical grid 
for the phase angle. In this case, the harmonic content of the solution is limited to 
one frequency with a wave length comparable to the computational domain which 
is resolved by 61 grid nodes. Therefore, all numerical schemes of concern provide 
accurate solutions. The numerical errors of |A| and 6 are tabulated in Table 2. 
There is slight advantage in using the higher-order schemes for the prediction of 
|A|; however, no obvious advantage of using the higher-order scheme for the phase 
angle calculation is observed. 

2.4.2 Shocked Sound Waves 

The second case is the propagation of shocked sound waves in a tube with a 
periodic boundary condition. The governing equations are the same as in the first 
case, namely, Eqns. (36) and (37), with cross section area (a.) equal to a constant. 
This case is interesting for its complex harmonic content compared to the first 
case. In addition, the capability of the high-order compact difference schemes for 
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shock capturing can also be studied. At time equal to zero, a sinusoidal pressure 
distribution is given. Because of the periodic boundary condition, only one cycle 
resolved by 61 grid nodes is imposed in the computational domain. According to 
the isentropic condition, the distributions of temperature, density, and speed of 
sound are also determined. The velocity profile is determined by the simple wave 
correlation 16 , such as 



dp 

p(x)c(x)' 


27 

7-1 



1 


(48) 


where the average flow properties are denoted by a bar. With the simple wave 
correlation, the wave forms of all flow properties are in phase and the initial con- 
dition of the present CFD computation matches the theoretical analysis provided 
by Morse and Ingard 8 . It is interesting to note that the simple wave correlation is 
an extension of a linear, plane, acoustic wave. For a variation of pressure less than 
5%, the plane wave relations could be adopted, such as 


r(l)= T(l + X_t^). 

M = <1 + ^ ) • <«> 

«<«) = ««)¥• 

where p'{x) is the prescribed pressure fluctuation. As shown in Eqns. (48) and (49), 
the wave speeds u + c,u-c and u vary as a result of the flow property distribution. 
The distortion of the wave form is a cumulative effect resulting from the wave speed 
distribution. For simple waves, i.e., all flow properties are in phase, the wave crest 
will quickly overtake the trough and form a shock. 

Figure 6 shows the time history of the pressure fluctuation at one end of the 
computational domain for various finite difference schemes. According to Morse and 
Ingard, the first shock appears after about two cycles for the case of a 10% pressure 
perturbation (p'/p = 0.1) 8 . All schemes of concern predict the wave steepening rate 
correctly. After the wave shocked, the flow evolution is no longer isentropic and the 
kinetic energy is gradually converted to thermal energy due to the existence of the 
shock wave. As a result, the strength of the shock wave diminishes as time passes. 

The shock front is a combination of many wave modes travelling at the same 
speed. The dispersion error introduced by the finite difference schemes will cause 
the high-wave-number waves to travel with erroneous speeds. As shown in Fig. 6, 
the methods of RK4-CD6 and RK4-CD4 with a small amount of the sixth-order 
artificial damping (77 = 0.02) crisply resolve the shock except for the over-shoots. 
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These over-shoots are caused by the Gibbs phenomenon and can be fixed only by 
TVD type shock-capturing schemes. Almost no difference can be observed between 
the results of the CD4 and CD6 methods. On the other hand, the method of RK4- 
CD6 without background filtering shows that significant high-wave-number waves 
lag behind the shock front because the compact difference scheme introduces no 
dissipative but high dispersive effects on the highest wave number waves. As shown 
in the figure, these oscillations eventually contaminate the whole solution. For the 
conventional RK4-CD2 method, results show significant oscillations of moderate 
wave numbers behind the shock front because of dispersion errors. 

Figure 7 shows the normalized power spectrums of the pressure profiles after 
about 17 cycles calculated by different methods. The analytical solution is plotted 
as the solid line. The power of each wave mode is roughly inversely proportional 
to the square of the wave number (oc 1/n 2 ). Since 61 grid nodes are used, only 
30 Fourier modes are resolved for the power spectrum (the other 30 modes are the 
complex conjugates). Clearly, the method of RK4-CD2 has significant errors in low 
wave modes. On the other hand, the methods of RK4-CD4 and RK4-CD6 compare 
well with the analytical solution. 

2.4.3 Vortex Propagation in an Uniform Flow 

A Lamb vortex propagated in an uniform flow is chosen as a two-dimensional 
numerical example. The vortex can be characterized by the circulation T and the 
core radius a. The azimuthal velocity uo at a distance r from the vortex center is 

given as, 


The flow near the vortex center is a rigid-body rotation ( u$ oc r). The flow far 
outside the core is irrotational (u e oc 1/r) with u e decreasing as r increases. Eqn. 
(50) is a continuous function to connect the two extremes. With the prescribed 
velocity field, the pressure and density distributions of the vortex can be determined 
by the momentum and the energy equations, 


dp _ tt| 
dr P r ’ 


7 V 




7~1 P 


(51) 

(52) 


where h a is the total enthalpy and is set to be a constant such as h c — 7P/(7 1)P 

with the free stream condition denoted by a bar. To proceed, substitute Eqns. (50) 
and (52) into Eqn. (51) and integrate the equation over r. As a result, the pressure 
distribution is obtained. Consequently, the density distribution and the whole flow 
field is determined. The solutions of this stationary vortex can be superimposed 
to any uniform flow with arbitrary speed. Physically, this process may be inter- 
preted as a stationary vortex being observed from a moving coordinate system with 
constant velocity. Thus, the vortex in an uniform flow can be constructed as, 
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u = u + u', 

v = V + v\ 

where the velocities of the background flow are denoted by a bar and the superscript 
' denotes the vortex velocities specified by Eqn. (50). The pressure and density 
distribution of the moving vortex is the same as that of the stationary vortex and 
may be obtained from the solutions of Eqns. (51) and (52). 

The boundary condition of the present case is an extension of the characteristic 
type treatment discussed in Case 1. Essentially, only one-dimensional character- 
istics (derived from two-dimensional flow equations) normal to the computational 
boundary are considered. For the purposes of this discussion, the subsonic out-flow 
condition is considered. The coupled equations of three out-running characteris- 
tics and one specified boundary condition, similar to that in Eqn. (47), should be 
solved. For steady state calculations, a back pressure (pb) is specified to regulate 
the flow rate, such as 


(53) 


k = 



(54) 


Similarly, the dimension of vectors Q and R is four and the matrix M -1 is a 4 x 4 
eigenvector matrix for the flux vector normal to the computational boundary. 


For a non-reflective boundary condition, Giles’ formulation 17 instead of the back 
pressure is used to fill the entry for the specified boundary condition, such as 


dc 4 
~dt 


+ (0,u,0, v) 


dy 



= 0 


(55) 


where y is in the direction parallel to the computational boundary. The variables 
c», i = 1, • • • , 4 are the characteristic variables and can be obtained by the similarity 
transformation from the non-conservative form equations as illustrated by Giles. In 
our case, the characteristic variables are derived from the conservative- form equa- 
tions using the same eigenvector matrix M -1 as in the aforementioned discussion. 
Giles’ non-reflective formulation, Eqn. (55), is relatively simple to used with an 
existing one-dimensional characteristic boundary condition. Nonetheless, according 
to Giles’ analysis, some two-dimensional effect is considered in the equation. Nu- 
merically, Eqn. (55) may be discretized according to the finite difference scheme 
of the interior nodes and combined with the discretized out-running characteristic 
equations (the two-dimensional version of Eqn. (46)) to form the complete subsonic 
out-flow boundary condition. It is interesting to note, however, according to Huff’s 
study 18 and our experience, that stretching the numerical grid nodes downstream 
to dampen the outgoing unsteady waves is just as effective. 
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For the in-flow conditions, the characteristic-type treatment combined with Giles 
equation (different from Eqn. (55)) may be adopted. For the present calculations, 
however, the upstream condition is relatively insensitive to various forms of non- 
reflective treatment as long as the proper out-running characteristic equation is 
selected and solved with the prescribed incoming conditions similar to that in Eqn. 
(47). In the present case, constant total pressure and total temperature are pre- 
scribed as the forcing boundary conditions upstream. 

As in the one-dimensional cases, dissipative and dispersive effects of various 
schemes are assessed. The prescribed vortex flow field contains a broad band of 
frequencies due to the distribution of the azimuthal velocity. Theoretically, all wave 
modes travels at the same speed to ensure the integrity of the vortex structure. 
For numerical methods with dispersive error, the shape of the vortex could deform, 
even break up in the later stage of the time marching procedure. In addition, the 
dissipation effect of finite difference schemes can be evaluated by the conservation 
of the sharp pressure dip at the center of the vortex propagated in the numerical 

grid. 

In the present calculations, the Mach number of the background flow is 0.4. The 
grid size is 301 x 91 in the streamwise and transverse directions, respectively. Uni- 
form grids are used in the axial direction and the transverse grids are stretched near 
the outside boundary. The CFL number calculated based on the background flow 
is about 0.7 for all calculations. As discussed before, a small amount of background 
filtering ( 77 = 0.02) is applied for all calculations. The core radius {a) is about 1 cm 
and is resolved by about 4 grid nodes. 

Figure 8 shows the vorticity and Mach number contours of the initial condition 
prescribed by Eqns. (50) - (53). Figure 9 shows the contours after the eddy prop- 
agates about 60 core radii downstream as simulated by various numerical schemes. 
The comparison between Fig. 8 and Fig. 9 shows that the structure of the eddy is 
retained by the compact difference schemes (CD4 and CD6). In contrast, the eddy 
predicted by the second-order central difference is shattered due to the excessive 

dispersive error. 

Figure 10 shows pressure distributions of the eddy at various instances. In this 
figure, the x axis represents the streamwise locations non-dimensionalized by the 
core radius of the vortex and the y axis is the pressure. For both CD4 and CD6 
methods, the pressure at the vortex center increases about 1 % through the process. 
In comparison, the results of the second-order scheme show pressure fluctuations 
with an overall increase about 3 %. The pressure fluctuation predicted by the 
second-order scheme is due to the deviation of the vortex path. 

2.4.4 Vortex Pairing 

Finally, the calculation of the single vortex is extended to the simulation of the 
vortex pairing. The vortex pairing is the controlling mechanism for the growth of 
the mixing layer 10 . In theory, vortex pairing occurs when the distance between 
two vortices is less than a threshold value. Unfortunately, no theoretical analysis 
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is available for compressible flows. In the present paper, the RK4-CD6 method is 
used to simulate the pairing process to demonstrate the resolution of the high-order 
compact difference scheme. 

The initial condition is specified by two identical vortices placed 5 core radii apart 
in a quiescent gas. The core radius is 1 cm and circulation is 15 m 2 /s. At the center 
of each vortex, there is a pressure deficit about 15% compared to the ambient gas. 
The grid size is 201 x 201. Uniform grids are used at the center of the computational 
domain to resolve the vortices. The grids are slightly stretched in all four directions 
to prevent erroneous wave reflection. In addition, one-dimensional characteristic 
equations combined with Giles’ unsteady, subsonic, out-flow equation is solved on 
the boundary as the non-reflective boundary condition. 

Figure 11 shows the contours of the vorticity magnitude at various stages of the 
vortex interaction. The whole sequence is about one and a half revolutions. Finally, 
a single larger vortex emerges as the result of the vortex pairing interaction. Figure 
12 shows the corresponding Mach number contours for the same flow. 

3. Concluding Remarks 

In this work, the quasi-one-dimensional and two-dimensional Euler solvers us- 
ing various combinations of the Runge Kutta methods and the compact difference 
schemes were developed for numerical simulations of unsteady flows. The accuracy 
of the finite difference schemes is assessed by Fourier analysis and numerical exam- 
ples in terms of numerical dissipation and dispersion. The dispersive characteristic 
is improved by high-order compact difference schemes compared to the second-order 
central difference. The increase of the order of time stepping scheme, on the other 
hand, enlarges the CFL limit for stable computations. In particular, significant im- 
provement of the dispersive effect is obtained by adopting the fourth-order compact 
scheme (6 to 8 grid nodes to resolve a wave) instead of the conventional second-order 
central difference (12 to 16 grid nodes for one wave). The use of the sixth-order 
compact scheme (5 to 8 grid nodes for one wave), however, gains little improvement 
compared with the fourth-order scheme. It was also found that the compact dif- 
ference schemes have no dissipative but high dispersive effects to the highest-wave- 
number waves resolved by a given numerical grid. Consequently, a small amount of 
background filtering is necessary to dissipate the high-wave-number waves and, at 
the same time, keep the low-wave-number solution intact. Other issues such as the 
order of accuracy of the Runge-Kutta schemes for nonlinear equations are analyzed. 
Specifically, the criteria for the 3 and 4-step methods to be third and fourth-order 
accurate are derived. The accuracy of the compact difference schemes for the steady 
state solution is also addressed. 

In general, simulation of unsteady flow provides an overwhelming amount of in- 
formation. It is our experience that the initial and boundary conditions must be 
carefully set up to obtain interpretable and physically meaningful solutions. For 
practical purposes, Giles’ non-reflective equations combined with one-dimensional 
characteristic equations and their implementation to the present numerical scheme 
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were illustrated in detail. In addition, the initial conditions of the simple wave, 
plane acoustic wave, and the Lamb vortex were also provided. Finally, 35 illus- 
trated in the numerical examples, for flows of simple harmonic content, e.g., one 
frequency in Case 1, the conventional second-order central difference scheme is ad- 
equate provided enough grid nodes are used to resolve the very wave mode. On 
the other hand, for flows of complex harmonic content, the use of the Runge-Kutta 
method combined high-order compact difference schemes shows crisp resolution of 
unsteady flows. 
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Table 2. The Relative Error of the Acoustic Admittance Calculation (%). 


Numerical Schemes Error of |A| Error of 6 
RK4-CD6 0.45 3.6 

RK4-CD4 0.52 3.3 

RK4-CD2 1.65 4.1 
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c. Dissipation of RK3-CD4 




e. Dissipation of RK3-CD2 



Fig. 1 Dissipation and dispersion characteristics of the RK3 time-stepping combined with 
various spatial discretization schemes. 
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CFL = 0.4 
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Fig. 2 Dissipation and dispersion characteristics of the RK4 time-stepping combined with 
various spatial discretization schemes. 
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Fig. 5 Acoustic admittance calculation using the RK4-CD6 method for the inlet perturbation 
as = 6. p'/p = 0.011. The solid line is Tsien’s theorem and the triangles are the 

CFD results, (a) Magnitude. (b)Phase angle. 
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Fig. 6 Time histories of the pressure fluctuations of the N-wave calculation at one end of the 
periodic domain by various numerical schemes. 
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(a) RK4— CD6 with A.D. wave mode 




(c) RK4— CD4 with A.D. wave mode 



(d) RK4-CD2 with A.D. wave mode 


Fig. 7 Power spectr um of the pressure distribution after about 17 cycles of the N-wa\e prop 
agation by various numerical schemes. 
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Fig. 10 Vortex pressure profiles at the center line at various instances predicted by different 
numerical schema 
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Fig- / / The vorticity magnitude contours for the vortex pairing. 
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Appendix B 

CMOTT Biweekly Seminars / Technical Meetings 


The purpose of these seminars is to exchange ideas and opinions on the latest 
developments and current state of turbulence and transition research. The speakers 
are invited from within and without of the NASA LeRC, including foreign speakers. 
The seminars were intended not noly to keep the members informed of the latest 
development of local turbulence and transition modeling research but also to increase 
interactions between group members and other researchers at the NASA LeRc. 

The following is the meeting schedule and the abstract of the semainars during 
the reporting period. 
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Date: July 3, 1991 

To: CMOTT Members and SVR and IFMD Staff 

From: William W. Liou (6682) 

Subject: CMOTT Biweekly Meeting 


The following is a tentative schedule for the CMOTT biweekly get-together from July 
10, 1991 to August 28, 1991. 

The presentations will be informal and active participation is expected from the at- 
tendants. Soda and snack will be served in the meetings. These meetings complement the 
CMOTT Seminar Series, which are mainly formal presentations. 

We would also appreciate some contributions from you. Subjects related to either the 
theoretical, experimental or computational aspects of turbulence and transition modeling 
are welcomed. Those who are willing to share their experience in these areas can contact 
me or Dr. T.-H. Shih at 6680 for further arrangement. 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 


July 10, 1991 

J. Lepicovsky (61-6753) 

LDV Measurement of Large Structures in a Tone 
Excited Turbulent Jet 

July 24, 1991 

C. R. Wang (5865) 

Computations of Turbulence in a Shock/Turbulent 
Boundary Layer Interaction Flow 

August 7, 1991 

A. Hsu (61-6648) 

PDF Turbulence Model and Its Applications 

August 28, 1991 

C. Steffen (8508) 

DTNS: An Accurate and Efficient Testbed for 
Incompressible Flow Turbulence Modeling 
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September 4, 1991 

CMOTT Members and SVR and IFMD Staff 
William W. Liou (6682) 

CMOTT Biweekly Meeting 

The following is a tentative schedule for the CMOTT biweekly get-together from 
September 11 , 1991 to October 23 , 1991 . 

The presentations will be informal and active participation is expected from the at- 
tendants. Soda and snack will be served in the meetings. These meetings complement the 
CMOTT Seminar Series, which are mainly formal presentations. 

We would also appreciate some contributions from you. Subjects related to either the 
theoretical, experimental or computational aspects of turbulence and transition modeling 
are welcomed. Those who are willing to share their experience in these areas can contact 
me or Dr. T.-H. Shih at 6680 for further arrangement. 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 



Date: 

To: 

From: 

Subject: 


Sept. 11, 1991 

T. Bui (5639) 

Implementation of the Chien Low-Re k-e Models into the 


Proteus Code 

Sept. 25, 1991 

K. Ahn (5965) 

A 2-D Oscillating Flow Analysis Using Quasi-steady 


Turbulence Model 

October 9, 1991 

J. Schwab (8446) 

Variable-density Turbulence Modeling for Turbomachinery 

October 23, 1991 

M. Mawid (5965) 

Multiphase Turbulent Combustion 
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Date: 

November 4, 1991 

To: 

CMOTT Members and SVR and IFMD Staff 

From; 

William W. Liou (3-6682) 

Subject: 

CMOTT Biweekly Meeting 


The following is a tentative schedule for the CMOTT biweekly get-together from 
November 6, 1991 to December 18, 1991. 

This will be the last session of the CMOTT group- meet ings/informal- seminars this 
year but the series will resume in mid-January 1992. Thank you for your patience and 
participation through out the year. The group-meetings/informal-seminars of CMOTT 
are meant not only to serve CMOTT members but also to provide an informal forum for 
those who are involved in transition/turbulence predictions. I thank all the speakers and 
participants who have made these objectives “realizable”. Now, we are planing for 1992. 
If you have any suggestions or like to give a talk or two in the coming year, please call me 
or Dr. T. H. Shih at 3-6680. In the mean time, don’t forget to mark your calendar for the 
following talks ! 


HAPPY HOLIDAYS Hi 

The meeting will start at 4:00 p.m. in Room 228, Sverdrup Building. 


Nov. 6, 1991 W. Liou (3-66S2) 

Weakly Nonlinear Models for Turbulent Free Shear Flows (2) 
- A Self-Contained Energy Transfer Model 

Nov 20, 1991 D. Ashpis (3-8317) 

DNS of Disturbances in Boundary Layer Flow 

Dec. 4, 1991 B. Rubinstein (61-6612) 

Analytical Theory of Turbulence and Turbulence Modeling: 
TSDIA and RNG 

Dec. 18, 1991 B. Duncan (61-2998) 

A New Three-Equation Model for Turbulence 
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Date: January 30, 1992 

To: CMOTT Members and SVR and IFMD Staff 


From: William W. Liou (3-6682) 

Subject: CMOTT Biweekly Meeting 


The following is a tentative program for the CMOTT biweekly get-together/seminar from 
February 5, 1992 to March 18, 1992. Ybu are welcomed to join us. 

Thanks to the your suggestions, we have made a few changes from last year’s format. First, 
we have scheduled the CMOTT Seminar Series, which are mainly formal presentations, into 
the biweekly time frame of the CMOTT group-meeting/informal-talks. Also, the abstract 
of each presentation, formal or informal, will be distributed about one week prior to its 
scheduled date. Again, if you are interested in giving a presentation, please contact us. 

The meeting will run from 1:30-2:15 PM in Room 145, Sverdrup Building, unless otherwise 
noted. 


(1) Feb. 5, 1992 D. Davis 

Weakly Nonlinear Vortex/Wave Interactions in 
Incompressible Cross-flow Boundary Layers in Transition 

(2) Feb. 19, 1992 Z. Yang 

A Modeling of Bypass Transition 

(3) March 4, 1992 K. Zaman 

Effect of ’’Delta- Tabs” on the Evolution of Axisymmetric 
Jets 

(4) March 18, 1992 Professor R. M. C. So, Arizona State University 

Near Wall Heat Transfer Modeling 
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Date: March 26, 1992 

To: CMOTT Members and SVR and IFMD Staff 

Prom: William W. Liou (3-6682) 

Subject: CMOTT Biweekly Meeting 

The following is a tentative program for the CMOTT biweekly get-together/seminar from 
April 1, 1992 to May 13, 1992. You are welcomed to join us. Also, if you are interested 
in giving a presentation, please let us know. 

The meeting will run from 1:30-2:15 PM in Room 228, Sverdrup Building, unless otherwise 
noted. 

(5) April 1, 1992 J. Van der Vegt 

The Development of an ENO-Osher Scheme for Direct 
Simulation of Compressible Flows 

(6) April 15, 1992 J. Goodrich 

Unsteady Time Asymptotic State: Incompressible Results, 

New Directions for Algorithms 

(7) April 29, 1992 T.-H. Shih 

Kolmogorov Behavior of Near- Wall Turbulence and 
Its Application in Turbulence Modeling 

(8) May 13, 1992 Z. Yang 

A Modeling of Bypass Transition 
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Date: June 1, 1992 

To: CMOTT Members and SVR and IFMD Staff 

From: William W. Liou (3-6682) 

Subject: CMOTT Biweekly Meeting 

The following is a tentative program for the CMOTT biweekly get-together/.seminar from 
June 10, 1992 to July 22, 1992. You are welcomed to join us. 

The talks will be informal and active participation will be expected from the audience. 

Also, if you are interested in giving a presentation about the progress or some results of 
your own work on turbulence or transition, please let us know. 

The meeting will run from 1:30-2:15 PM in Room 228, Sverdrup Building, unless otherwise 
noted. 

(9) June 10, 1992 J. Zhu 

Finite Volume Computations in Incompressible Flows 
with Complex Geometries 

(10) June 24, 1992 J. Lee 

RPLUS Code and Standard k - c Models and Applications 

(11) July 8, 1992 R. Mankbadi 

Unsteady Turbulent Flows 

(12) July 22, 1992 D. Rigby 

The Effect of Spanwise Variations in Momentum on 
Leading Edge Heat Transfer 
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Biweekly Meeting Series (1) 


Weakly Nonlinear Vortex/Wave Interactions in 
Incompressible Crossflow Boundary Layers in 

Transition 

by 

Dominic Davis 

ICOMP 


Wed., 5 February, 1992 
1:30-2:15 PM 
Room 145, SVR Building 


ABSTRACT 

The instability of an incompressible three-dimensional boundary layer 
is considered theoretically and computationally in the context of vor- 
tex/wave interactions. Specifically the work centres on two low-amplitude, 
lower-branch Tollmien-Schlichting (TS) waves which mutually interact 
to induce a weak longitudinal vortex flow;the vortex motion, in turn.gives 
rise to significant wave-modulation via wall-shear forcing.The character- 
istic Reynolds number is taken as a large parameter and.as a conse- 
quence.the TS waves and the vortex are governed primarily by triple- 
deck theory.The nonlinear interaction is captured by a viscous partial- 
differential system for the vortex coupled with a pair of amplitude equa- 
tions for the wave pressures. Computations were performed for relatively 
small crossflow values.Three distinct possibilities were found to emerge 
for the nonlinear behaviour of the flow solution downstream - an alge- 
braic finite-distance singularity.far- downstream decay or repeated oscil- 
lations - depending on the various parameter values.the input amplitudes 
and the wave angles. 
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Biweekly Meeting Series (1992-2) 

A k-e Model for Near Wall Turbulence and 
its Application in Turbulent Boundary Layer 
with/without Pressure Gradient 

by 

Z. Yang 

ICOMP 


Wed., 19 February, 1992 
1:30-2:15 PM 
Room 145, SVR Building 


ABSTRACT 

A fc - c model is proposed for turbulent wall bounded flows. In this 
model, turbulent velocity scale and turbulent time scale are used to 
define the eddy viscosity. The time scale used is bounded from below 
by the Kolmogorov time scale. The dissipation equation is reformulated 
using this time scale, removing the singularity of the high Reynolds 
number J fc - t equation at the wall and rendering the introduction of 
the pseudo-dissipation unnecessary. The damping function used in the 
eddy viscosity is chosen to be a function of Ry = instead of y + . 
Thus, the model could be used for flows with separation. The model 
constants used are the same as the model constants in the commonly 
used high turbulent Reynolds standard k-e model. Thus, the proposed 
model would reduce to the standard k-e model when it is far away 
from the wall. Boundary layer flows at zero pressure gradient, favorable 
pressure gradient, adverse pressure gradient and increasingly adverse 
pressure gradient are calculated respectively. The comparisons of model 
predictions and the available experimental data are found to be good. 
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Biweekly Meeting Series (1992-3) 


Effect of Tabs on the Evolution of Axisymmetric Jets 

by 

Khairul Zaman 

IFMD 

Wed., 4 March, 1992 
1:30-2:15 PM 
Room 145, SVR Building 


ABSTRACT 

Vortex generators, in the form of small tabs at the nozzle exit, can have 
a profound influence on the evolution of an axisymmetric jet. Using 
tabs of certain shapes, the jet cross section can be distorted almost 
arbitrarily. Such distortion is accompanied by elimination of screech 
noise from supersonic jets and a significant increase in jet mixing. Key 
results obtained so far, covering a jet Mach number range of 0.3 and 
1.8, will be summarized in this presentation. Observations will be made 
on the mechanisms of the effect including the likely vorticity dynamics 
in the flow. 
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Biweekly Meeting Series (1992-4) 


Near-Wall Modeling of Turbulent Heat Transfer 

by 

Professor Ronald M. C. So 

Mechanical and Aerospace Engineering 
Arizona State University 


Wed., 18 March, 1992 
1:30-2:30 PM 
Room '119, Building 5 


ABSTRACT 

A near-wall two-equation model for turbulent heat fluxes is derived from 
the temperature variance and its dissipation-rate equations and the as- 
sumption of gradient transport. The near-wall asymptotics of each term 
in the exact equations are examined and used to derive near-wall cor- 
rection functions that render the modeled equations consistent with 
these behavior. Thus modeled, the equations are used to calculate 
fully-developed pipe and channel flows with hear transfer. It is found 
that the proposed two-equation model yields asymptotically correct near- 
wall behavior for the normal heat flux, the temperature variance and its 
near-wall budget and correct limiting wall values for these properties 
compared to direct simulation data and measurements obtained under 
different wall boundary conditions. 


CONTACT: T.-H. Shih, PABX 3-5698 
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Biweekly Meeting Series (1992-5) 


The Development of an ENO-Osher Scheme for 
Direct Simulation of Compressible Flows 

by 

Jaap Van der Vegt 
ICOMP 


Wed., April 1, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 


Direct simulation of turbulence and transition in compressible wall 
bounded flows presents an alternative to investigate important physical 
phenomena which are difficult to measure or study otherwise. It also 
provides data useful for turbulence modeling. A 1 new program is be- 
ing developed which solves the three-dimensional compressible Navier- 
Stokes equations using a higher order, fully implicit and time accurate 
ENO scheme together with Osher flux splitting. In this presentation an 
overview will be given of the numerical scheme and several test cases, 
both for supersonic and subsonic flow, will be presented and further 
improvements will be discussed. 
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Biweekly Meeting Series (1992-6) 

Unsteady Time Asymptotic States: Incompressible Results 
and New Directions for Algorithms 

by 

John Goodrich 

IFMD 

Wed., April 15, 1992 
1:30-2:15 PM 
Room 228, SVR Building 

ABSTRACT 

Unsteady time asymptotic flow states for high Reynolds number viscous incom- 
pressible flow problems are presented. Discrete frequency flows are shown for 
the square driven cavity, with periodic cases for Re = 9000 and Re - 9600 , and 
with aperiodic cases for Re = 9700 and Re = 10000 . The algorithm for these cal- 
culations is based on the fourth order PDE for incompressible fluid flow which 
uses the streamfunction as the only dependent variable, it is second order ac- 
curate in space and time, and it has a stability constraint of CFL < l. I he 
algorithm is extremely robust with respect to Reynolds number. 

The direct numerical simulation of tr ansition and turbulence requires nu- 
merical methods to be more than second order accurate in order to accurately 
represent the relevant scales of the physical processes. Recently deveope 
finite difference algorithms are presented for unsteady convection equations, 
including the advection and inviscid Burgers equation in one space dimension, 
and the wave equation treated as a system, with remarks on diffusion equa- 
tions and extension to higher space dimensions. The new algorithms that will 
be discussed all use local stencils, they range from third to seventh order in 
accuracy, they all have the same order of accuracy in both space and time, and 
they are all one step explicit methods (except for diffusion equations). Since all 
of the algorithms use a small local stencil, the number of degrees of freedom of 
known data required for higher order accuracy is obtained by higher information 
density than just the solution data. The use of a two point stencil (for sonie 
of the methods) allows for arbitrary grid spacing, though a convective stability 
constraint must be observed at each grid point. The use of local data for an 
explicit algorithm with high order accuracy avoids the requirement of using a 
global solution method such as compact differencing or spline based algorithms. 
There will be computational results for all of the algorithms that are presented. 




pjlCOMPr 
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Biweekly Meeting Series (1992-7) 

Kolmogorov Behavior of Near-Wall Turbulence 
and Its Application in Turbulence Modeling 

by 

Tsan-Hsing Shih 

ICOMP 



Wed., April 29, 1992 
1:30-2:15 PM 
Room 228, SVR Building 


ABSTRACT 

The near-wall behavior of turbulence is re-examined in a way different from 
that proposed by Hanjalic and Launder^ and followersI 2 W 3 l’l 4 ll 5 '. It is shown 
that at a certain distance from the wall, all energetic large eddies will reduce 
to Kolmogorov eddies (the smallest eddies in turbulence). All the important 
wall parameters, such as friction velocity, viscous length scale, and mean strain 
rate at the wall, are characterized by Kolmogorov microscales. According to this 
Kolmogorov behavior of near-wall turbulence, the turbulence quantities, such 
as turbulent kinetic energy, dissipation rate, etc. at the location where the 
large eddies become "Kolmogorov" eddies, can be estimated by using both di- 
rect numerical simulation (DNS) data and asymptotic analysis of near-wall 
turbulence. This information will provide useful boundary conditions for the 
turbulent transport equations. As an example, the concept is incorporated in 
the standard k-e model which is then applied to channel and boundary layer 
flows. Using appropriate boundary conditions (based on Kolmogorov behavior 
of near-wall turbulence), there is no need for any wall-modification to the k- 
e equations (including model constants). Results compare very well with the 
DNS and experimental data. 
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Biweekly Meeting Series (1992-8) 


A Modeling of Bypass Transition 


by 

Z. Yang 

ICOMP 


Wed., May 13, 1992 
1:30-2:15 PM 
Room 228, SVR Building 


ABSTRACT 

A model for the calculation of bypass transitional boundary layers due to the 
freestream turbulence is proposed. The model combines a near wall k-e model 
proposed for the fully developed turbulent flows with the intermittency of the 
transitional boundary layers. The intermittency factor is assumed to be a func- 
tion of both the free stream turbulence and the shape factor of the boundary 
layer. Transitional boundary layers over a flat plate with different freestream 
turbulence level are calculated using the proposed model. It is found that the 
model calculations agree well with the experimental data and give a better pre- 
diction compared with other low Reynolds number k - e models, which do not 
incorporate the intermittency effect. 
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Biweekly Meeting Series (1992-9) 

Finite-Volume Computations of Incompressible 
Flows with Complex Geometries 

by 



J. Zhu 

ICOMP 


Wed., June 10, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 

A brief review is given of finite-volume procedures developed at the 
Institute for Hydromechanics, University of Karlsruhe, for calculating 
incompressible elliptic flows with complex boundaries. The procedures 
include: numerical grid generation, higher-order bounded convection 
schemes, zonal solution, simulation of two-phase flows, and near-wall 
turbulence modelling. Various application examples will be given. 
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Biweekly Meeting Series (1992-10) 

Development of the RPLUS code with 
Standard *-* model and Their Applications 

by 

Jinho Lee 

Sverdrup Technology, Inc. 


Wed., June 24, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 


The primary goal of this research effort is to develop a CFD tool which can be used in 
a variety of practical supersonic/hypersonic propulsion device development /analysis envi- 
ronments, One focus of this work has been to develop and validate the Reactive Propulsion 
code based on LU Scheme(RPLUS). This effort also includes the development of turbu- 
lence models which can be used in the predictions of highly complex flow environments 
inside of combustors. 

This presentation will cover only a small part of a larger development effort and focus 
will primarily on the analysis, implementation, and development of the turbulence model 
capabilities of the RPLUS code. 

Some of the issues which will be covered are; formulation of the turbulence models, 
the numerical technique used to solve the turbulence model equations, and modeling of 
compressibility effects. The primary numerical technique used in the RPLUS code is 
the LU-SSOR(LU scheme based on Successive Symmetric Over Relaxation) technique. 
Therefore, the turbulent kinetic energy transport and dissipation transport equations are 
also solved using the LU-SSOR numerical technique. 

Both two and three dimensional turbulence model development are being developed 
for the RPLUS code. However, the majority of the presentation will focus on the devel- 
opment of the two dimensional k-e models for the RPLUS code. Issues regarding com- 
pressible wall- function boundary conditions and compressibility effects will be addressed. 
Both low and high Reynolds number forms of the k-c models are being developed. The 
“standard” low Reynolds number model of Launder- Sharma ?md Chien has been used in 
this study. The problems of primary interests are supersonic turbulent boundary- layers, 
shock- wave /bound ary- layer interactions, and shear-layers in two or three dimensional en- 
vironments. 
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Biweekly Meeting Series (1992-11) 


Unsteady Turbulent Flows 

by 

Reda R. Mankbadi 

NASA Lewis Research Center 


Wed., July 8, 1992 
1:30-2:30 PM 
Room 228, SVR Building 

ABSTRACT 

Current research activities emphasize computation/modelling of turbu- 
lent flows when basic flow is time-periodic. Wall-bounded flows and free 
shear flows exhibit different features when the basic flow is unsteady; 
and as such, different approaches are used to model them. 

(A) In wall-bounded, oscillatory flows, two approaches are used to 
model turbulence: (I) Turbulence is assumed to behave in a quasi-steady 
manner, and steady-state models are directly extended to the unsteady 
case. This approach fails at high frequencies of oscillations. (II) Rapid 
distortion theory (RDT) is successfully adapted to aid in turbulence 
modelling of highly unsteady flows (high frequencies). The eddy viscosity 
hypothesis is replaced by the ratio of turbulent stresses/kinetic energy; 
which is given by RDT as a function of the accumulated rate of strain. 

(Bl In free shear flows (naturally unsteady, or excited to be un- 
steady), two approaches are investigated: (I) The large-scale (organized, 
coherent) component is modelled as instability waves interacting with 
each other as well as with the mean flow and the fine-scale (random, 
background) turbulence. Integrated kinetic-energy equations are then 
obtained for each scale of motion. The approach is successful in pre- 
dicting results in good agreements with experiments in which excitation 
devices are used to control jet mixing and turbulence. (II) The other 
approach adopted is Large-Eddy Simulations (LES) with application to 
predicting the far-field noise of a supersonic jet. 


CONTACT: William W. Liou. PABX 3-C682 
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Biweekly Meeting Series (1992-12) 


The Effect of Spanwise Variations in 
Momentum on Leading Edge Heat Transfer 

by 

David Rigby 

Sverdrup Tech. INc. 


Wed., July 22, 1992 
1:30-2:30 PM 
Room 228, SVR Building 


ABSTRACT 

A study of the effect of spanwise variation in momentum on leading 
edge heat transfer is discussed. Numerical and experimental results 
are presented for a circular leading edge and for a 3:1 elliptical leading 
edge. Direct comparison of the two-dimensional results, that is with no 
spanwise variations, to the analytical results of Frossling is very good. 
The numerical calculation, using the PARC3D code, solves the three- 
dimensional Navier-Stokes equations, assuming steady laminar flow on 
the leading edge region. Experimentally, increases in spanwise averaged 
heat transfer coefficient as high as 50% above the two-dimensional value 
were observed. Numerically, the heat transfer coefficient was seen to 
increase by as much as 25% percent. In general, the circular leading 
edge under the same flow conditions, produced a higher heat transfer 
rate than the elliptical leading edge. As a percentage of the respective 
two-dimensional values, the circular and elliptical leading edges showed 
similar sensitivity to spanwise variations in momentum. By equating the 
root mean square of the amplitude of the spanwise variation in momen- 
tum to the turbulence intensity, a qualitative comparison between the 
present work and turbulent results was possible. 


CONTACT: William W. Liou, PABX 3-6682 
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